#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 17 23:13:46 2021

@author: liqingsimac
"""

'''
###########################################
#https://stackoverflow.com/questions/21352580/
#matplotlib-plotting-numerous-disconnected-line-segments
#-with-different-colors

import numpy as np
import pylab as pl
from matplotlib import collections  as mc

lines = [[(0, 1), (1, 1)], [(2, 3), (3, 3)], [(1, 2), (1, 3)]]
c = np.array([(1, 0, 0, 1), (0, 1, 0, 1), (0, 0, 1, 1)])

lc = mc.LineCollection(lines, colors=c, linewidths=2)
fig, ax = pl.subplots()
ax.add_collection(lc)
ax.autoscale()
ax.margins(0.1)
'''


###########################################
##Example 1-1-自由落体

import numpy as np
import matplotlib.pyplot as plt

g=9.8; C2=100; C1=[0,10,20,30,40]; t=np.linspace(0,12,101)

fig=plt.figure(); ax=fig.add_subplot(111)
for k in range(len(C1)):
    y=-g*t**2/2+C1[k]*t+C2
    ax.plot(t,y)

ax.set_xlim([0,12]); ax.set_ylim([0,240])
ax.set_xlabel('t=time'); ax.set_ylabel('y=height')
ax.set_title('Falling object with different initial velocity upwards')
ax.annotate(r'initial height $C_2=100$',xy=(0.5,220))
ax.annotate(r'initial upwards velocity $C_1=0,10,20,30,40$',xy=(0.5,200))
ax.annotate(r'$y=-\frac{1}{2}gt^2+C_1t+C_2$',xy=(7,160))
#fig.savefig('ode-example-1-1-0.png')


'''
###########################################
##Example 1-1-2
import numpy as np
import matplotlib.pyplot as plt

t=np.linspace(0,2*np.pi,31)
x0=np.arange(-3,3,step=1)
y0=np.arange(-3,3,step=1)

fig=plt.figure()
ax=fig.add_subplot(111)

for k in range(len(x0)):
    for m in range(len(y0)):
        flipcoin=np.random.rand()
        if flipcoin<0.2:
            r=np.sqrt(x0[k]**2+y0[m]**2)
            x=x0[k]+r*np.cos(t)
            y=y0[m]+r*np.sin(t)
            ax.plot(x,y)
        else:
            continue
ax.set_xlim([-10,10])
ax.set_ylim([-10,10])
ax.set_aspect('equal')
#fig.savefig('ode-example-1-1-2.png')
'''

'''
###########################################
##Example 1-1-2-second version
import numpy as np
import matplotlib.pyplot as plt

N=6 # number of circles
M=10 # maximal xy limit of the center of the circle
t=np.linspace(0,2*np.pi,31)

fig=plt.figure()
ax=fig.add_subplot(111)

for k in range(N):
    x0=M*np.random.rand()-5
    y0=M*np.random.rand()-5
    r=np.sqrt(x0**2+y0**2)
    x=x0+r*np.cos(t)
    y=y0+r*np.sin(t)
    ax.plot(x,y)
    
ax.set_xlim([-10,10])
ax.set_ylim([-10,10])
ax.set_aspect('equal')
ax.set_yticks([-10,-5,0,5,10])
#fig.savefig('ode-example-1-1-2.png')
'''

'''
###########################################
##Example 1-1-2-third version 很多圆
import numpy as np
import matplotlib.pyplot as plt

N=10 # number of circles
t=np.linspace(0,2*np.pi,31)

fig=plt.figure()
ax=fig.add_subplot(111)

for k in range(N):
    x0=20*np.random.rand()-10
    y0=20*np.random.rand()-10
    r=5*np.random.rand()
    x=x0+r*np.cos(t)
    y=y0+r*np.sin(t)
    ax.plot(x,y)
    
ax.set_xlim([-10,10])
ax.set_ylim([-10,10])
ax.set_aspect('equal')
ax.set_xticks([])
ax.set_yticks([])
'''

'''
###########################################
## 1-1-4 常微分方程的解
import numpy as np
import matplotlib.pyplot as plt

fig=plt.figure()
ax=fig.add_subplot(111)

x1=np.linspace(0.01,2.9,31)
y1=1/x1 + x1**4*(1/5)
x2=np.linspace(-2.9,-0.01,31)
y2=1/x2 + x2**4*(1/5)
ax.plot(x1,y1,'r-')
ax.plot(x2,y2,'b-')
ax.set_xlim([-3,3])
ax.set_ylim([-20,20])
#ax.set_aspect('equal')
#ax.set_xticks([])
ax.set_yticks([-20,-10,0,10,20])
fig.savefig('ode-1-1-4-solution.png')
'''

'''
###########################################
##例子1-1-Jacobi行列式 version 1

import numpy as np
import matplotlib.pyplot as plt

[X,Y] = np.mgrid[-2:2:17j,-2:2:17j]

mycolor=[['b','g','r'],['m','y','k']]
mylinestyle=[['-','--','-.'],['-','--','-.']]
fig=plt.figure()
ax1=fig.add_subplot(1,2,1)

# =============================================================================
for k in (0,3,6,9,12,15):
    ax1.plot(X[k,:],Y[k,:])
# for m in (0,1,2):
#     ax1.plot(X[:,0],Y[:,8*m],color=mycolor[1][m],linestyle=mylinestyle[1][m])
# ax1.set_xlim([-3,3])
# =============================================================================

#ax1.plot(X[0,:],Y[0,:],color=mycolor[0][0],linestyle=mylinestyle[0][0])
#ax1.plot(X[4,:],Y[4,:],color=mycolor[0][0],linestyle=mylinestyle[0][0])
#ax1.plot(X[8,:],Y[8,:],color=mycolor[0][1],linestyle=mylinestyle[0][1])
#ax1.plot(X[12,:],Y[12,:],color=mycolor[0][1],linestyle=mylinestyle[0][1])
#ax1.plot(X[16,:],Y[16,:],color=mycolor[0][2],linestyle=mylinestyle[0][2])
#ax1.plot(X[:,0],Y[:,0],color=mycolor[1][0],linestyle=mylinestyle[1][0])
#ax1.plot(X[:,8],Y[:,8],color=mycolor[1][1],linestyle=mylinestyle[1][1])
#ax1.plot(X[:,16],Y[:,16],color=mycolor[1][2],linestyle=mylinestyle[1][2])


ax1.set_xlim([-3,3])
ax1.set_ylim([-3,3])
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_aspect('equal')
ax1.set_xticks([-2,0,2])
ax1.set_yticks([-2,0,2])

U=X**2-Y**2
V=X*Y

ax2=fig.add_subplot(1,2,2)

# =============================================================================
for k in (0,3,6,9,12,15):
    ax2.plot(U[k,:],V[k,:])
# for m in (0,1,2):
#     ax2.plot(U[:,0],V[:,8*m],color=mycolor[1][m],linestyle=mylinestyle[1][m])
# ax2.set_xlim([-5,5])
# =============================================================================

#ax2.plot(U[0,:],V[0,:],color=mycolor[0][0],linestyle=mylinestyle[0][0])
#ax2.plot(U[4,:],V[4,:],color=mycolor[0][0],linestyle=mylinestyle[0][0])
#ax2.plot(U[8,:],V[8,:],color=mycolor[0][1],linestyle=mylinestyle[0][1])
#ax2.plot(U[12,:],V[12,:],color=mycolor[0][1],linestyle=mylinestyle[0][1])
#ax2.plot(U[16,:]+0.3,V[16,:],color=mycolor[0][2],linestyle=mylinestyle[0][2])
#ax2.plot(U[:,0],V[:,0],color=mycolor[1][0],linestyle=mylinestyle[1][0])
#ax2.plot(U[:,8],V[:,8],color=mycolor[1][1],linestyle=mylinestyle[1][1])
#ax2.plot(U[:,16]-0.3,V[:,16],color=mycolor[1][2],linestyle=mylinestyle[1][2])

ax2.set_xlim([-5,5])
ax2.set_ylim([-5,5])
ax2.set_xlabel(r'$u=x^2-y^2$')
ax2.set_ylabel(r'$v=xy$')
ax2.set_aspect('equal')
ax2.set_xticks([-4,0,4])
ax2.set_yticks([-4,0,4])

fig.tight_layout(w_pad=5)
#fig.savefig('ode-example-1-1-jacobi.png')
'''

'''
###########################################
##例子1-1-Jacobi行列式 version 2

import numpy as np
import matplotlib.pyplot as plt

[X,Y] = np.mgrid[-2:2:17j,-2:2:17j]
U=X**2-Y**2
V=X*Y

fig=plt.figure()

###########################################
ax1=fig.add_subplot(2,2,1)
for k in (0,3,6,9,12,15):
    ax1.plot(X[k,:],Y[k,:])

ax1.set_xlim([-3,3])
ax1.set_ylim([-3,3])
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_aspect('equal')
ax1.set_xticks([-2,0,2])
ax1.set_yticks([-2,0,2])

###########################################
ax2=fig.add_subplot(2,2,2)
for k in (0,3,6,9,12,15):
    ax2.plot(U[k,:],V[k,:])

ax2.set_xlim([-5,5])
ax2.set_ylim([-5,5])
ax2.set_xlabel(r'$u=x^2-y^2$')
ax2.set_ylabel(r'$v=xy$')
ax2.set_aspect('equal')
ax2.set_xticks([-4,0,4])
ax2.set_yticks([-4,0,4])

###########################################
ax3=fig.add_subplot(2,2,3)
for k in (0,3,6,9,12,15):
    ax3.plot(X[:,k],Y[:,k])

ax3.set_xlim([-3,3])
ax3.set_ylim([-3,3])
ax3.set_xlabel('x')
ax3.set_ylabel('y')
ax3.set_aspect('equal')
ax3.set_xticks([-2,0,2])
ax3.set_yticks([-2,0,2])

###########################################
ax4=fig.add_subplot(2,2,4)
for k in (0,3,6,9,12,15):
    ax4.plot(U[:,k],V[:,k])

ax4.set_xlim([-5,5])
ax4.set_ylim([-5,5])
ax4.set_xlabel(r'$u=x^2-y^2$')
ax4.set_ylabel(r'$v=xy$')
ax4.set_aspect('equal')
ax4.set_xticks([-4,0,4])
ax4.set_yticks([-4,0,4])

###########################################

fig.tight_layout(w_pad=-5,h_pad=2)
#fig.savefig('ode-example-1-1-jacobi.png')
'''

'''
###########################################
##例子1-1-Jacobi行列式 version 3

import numpy as np
import matplotlib.pyplot as plt

[X,Y] = np.mgrid[-2:2:17j,-2:2:17j]
U=X**2-Y**2
V=X*Y

mycolor=['b','g','r','m','y','k']

fig=plt.figure()

###########################################
ax1=fig.add_subplot(1,2,1)
for k in range(6):
    ax1.plot(X[3*k,:],Y[3*k,:],color=mycolor[k])
    ax1.plot(X[:,3*k],Y[:,3*k],color=mycolor[k])
    
ax1.set_xlim([-3,3])
ax1.set_ylim([-3,3])
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_aspect('equal')
ax1.set_xticks([-2,0,2])
ax1.set_yticks([-2,0,2])

###########################################
ax2=fig.add_subplot(1,2,2)
for k in range(6):
    ax2.plot(U[3*k,:],V[3*k,:],color=mycolor[k])
    ax2.plot(U[:,3*k],V[:,3*k],color=mycolor[k])

ax2.set_xlim([-5,5])
ax2.set_ylim([-5,5])
ax2.set_xlabel(r'$u=x^2-y^2$')
ax2.set_ylabel(r'$v=xy$')
ax2.set_aspect('equal')
ax2.set_xticks([-4,0,4])
ax2.set_yticks([-4,0,4])

fig.tight_layout(w_pad=3)
fig.savefig('ode-example-1-1-jacobi-2.png')
'''

'''
###########################################
## Example 1-2-0-a 积分曲线族 y=tan(x+C)
import numpy as np
import matplotlib.pyplot as plt

x=np.linspace(-np.pi/2+0.01, np.pi/2-0.01,31)
y=np.tan(x)
C=np.linspace(-3,3,7)

fig=plt.figure()
ax=fig.add_subplot(111)

for k in range(len(C)):
    ax.plot(x+C[k],y,'-')

ax.set_xlim([-5,5])
ax.set_ylim([-20,20])
ax.set_xlabel('x')
ax.set_ylabel('y')
fig.savefig('ode-example-1-2-0-a.png')
'''

'''
###########################################
##线素场 Example 1-2-0-b Riccati Equation 
import numpy as np
import matplotlib.pyplot as plt

def myfun(x,y):
    return x**2+y**2
    
stk=0.2 # length of each stick (line element)

x=np.linspace(0,2,num=5)
y=np.linspace(0,2,num=5)

fig=plt.figure()
ax=fig.add_subplot(111)
for k in range(5):
    for m in range(5):
        xp=x[k]
        yp=y[m]
        slope=myfun(xp,yp)
        dx=stk/np.sqrt(1+slope**2)
        dy=dx*slope
        ax.plot([xp-dx/2,xp+dx/2],[yp-dy/2,yp+dy/2],color='b')

ax.set_aspect('equal')
fig.savefig('ode-example-1-2-0-b.png')
'''


'''
###########################################
##线素场 Example 1-2-1 
import numpy as np
import matplotlib.pyplot as plt

def myfun(x,y):
#    if (y==0):
#        return np.inf
#    else:
#        return y/x
#        return -x/y
#    return x**2+y**2
    a=5
    numerator=y*((x-a)**2+y**2)**(3/2) - y*((x+a)**2+y**2)**(3/2)
    denominator=(x+a)*((x-a)**2+y**2)**(3/2) - (x-a)*((x+a)**2+y**2)**(3/2)
    return numerator/denominator
    
b=10; N=21; 
stk=0.3 # length of each stick (line element)

x=np.linspace(-b,b,num=N)
y=np.linspace(-b,b,num=N)

fig=plt.figure()
ax=fig.add_subplot(111)
for k in range(N):
    for m in range(N):
        xp=x[k]
        yp=y[m]
        slope=myfun(xp,yp)
        if (slope==np.inf):
            ax.plot([xp,xp],[yp-stk/2,yp+stk/2])
        else:
            dx=stk/np.sqrt(1+slope**2)
            dy=dx*slope
            ax.plot([xp-dx/2,xp+dx/2],[yp-dy/2,yp+dy/2])

ax.set_aspect('equal')
fig.savefig('ode-example-1-2-1.png')
'''
'''
###########################################
##线素场 期中考试#2 
import numpy as np
import matplotlib.pyplot as plt

def myfun(x,y):
    return x*y-1
    
stk=0.3 # length of each stick (line element)

x=np.array([0,1,2])
y=np.array([0,1,2])

fig=plt.figure()
ax=fig.add_subplot(111)
for k in range(len(x)):
    for m in range(len(y)):
        xp=x[k]
        yp=y[m]
        slope=myfun(xp,yp)
        dx=stk/np.sqrt(1+slope**2)
        dy=dx*slope
        ax.plot([xp-dx/2,xp+dx/2],[yp-dy/2,yp+dy/2],
                color='b',linewidth=3)

ax.set_aspect('equal')
ax.set_xlim([-1,3])
ax.set_ylim([-1,3])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.hlines(y=0,xmin=-1,xmax=3,colors='b')
ax.vlines(x=0,ymin=-1,ymax=3,colors='b')
ax.set_xticks([-1,0,1,2,3])
ax.set_yticks([-1,0,1,2,3])

fig.savefig('ode-midterm-problem-2.png')
'''





